First Principles Calculations of Shock Compressed Fluid Helium 
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The properties of hot dense helium at megabar pressures were studied with two first-principles 
computer simulation techniques, path integral Monte Carlo and density functional molecular dy- 
namics. The simulations predicted that the compressibility of helium is substantially increased by 
electronic excitations that are present in the hot fluid at thermodynamic equilibrium. A maximum 
compression ratio of 5.24(4)-fold the initial density was predicted for 360 GPa and 150 000 K. This 
result distinguishes helium from deuterium, for which simulations predicted a maximum compression 
ratio of 4.3(1). Hugoniot curves for statically precompressed samples are also discussed. 



There has been considerable controversy in the deu- 
terium equation of state (EOS) since laser shock wave 
experiments probed the megabar pressure regime for the 
first time and predicted that deuterium is highly com- 
pressible under shock conditions to approximately 6-fold 
the initial density 0, 0] ■ Such a high compression ratio 
was neither reproduced with magnetic compression ex- 
periments 0, Q nor with explosively driven shocks [jj . 
Both sets of later experiments predicted compression ra- 
tios close to 4.3(1), which is in good agreement with re- 
sults from first principles computer simulations 0, 0, . 
When we applied the same simulation techniques, path 
integral Monte Carlo (PIMC) and density functional 
molecular dynamics (DFT-MD), to hot dense helium, we 
found that helium's shock compressibility is substantially 
increased due to electronic excitations in the fluid. 

In this Letter, we make the prediction that electronic 
excitations in helium lead to a maximum shock compres- 
sion ratio of 5.24(4), while such excitations in deuterium 
do not increase the compressibility ratio beyond 4.3(1). 
Furthermore, we show that the compression ratio is re- 
duced when the sample is precompressed statically in a 
diamond anvil cell before a shock is launched. Such novel 
compression techniques are currently developed and data 
for dense helium are forthcoming [lOj. The combination 
of static and dynamic compression techniques allows the 
study of materials at much higher densities, and their 
application to hydrogen and helium will enable a direct 
characterization of a much larger section of the isentrope 
that determines the interiors of giant planets. 

Present studies of giant planetary interiors [Tlj are 
based on approximate free energy models that rely on 
analytical thermodynamic expressions and are often fit 
to experimental results if available. Although these mod- 
els (for helium read 0[H0) 

are very practical, their 
predictive capabilities for shock states and the EOS are 
limited because interaction and polarization effects in a 
dense fluid are very difficult to study analytically, which 
underlines the need of first principles simulations. 

Shock wave experiments provide us with direct infor- 
mation for materials' EOS at high pressure and tempera- 
ture. When a shock wave passes through the sample, the 



thermodynamic state of the material, characterized by 
the internal energy, pressure, and volume, changes from 
initially (E ,P , V ) to the final values of (E, P, V). The 
conservation of mass, momentum, and energy yields the 
Hugoniot condition [15J, 



H = (E -Eq) + -(P + P )(V - V ) = 0. 



(1) 



The resulting Hugoniot curve is the locus of all final 
states that can be reached for different shock velocities. 
Theoretically, the Hugoniot curve can be calculated from 
the EOS, assuming Pq <C P and using the initial molec- 
ular volume from the experiment, Vq — 32 A cm 3 /mol 
(po = 0.1235 gem -3 ) 0. For E one takes the energy 
of an isolated helium atom, which must be consistent 
with the internal energy E(V, T) derived from a partic- 
ular method. For PIMC we use E = -79.0048 eV per 
atom because PIMC is exact for the helium atom. The 
assumptions for Po and Eq also remain sufficiently valid 
for samples that have been precompressed statically to 4 
times po- 

PIMC is a finite-temperature quantum simulation 
method that we used to model dense fluid helium as a sys- 
tem of nuclei and electrons that interact via the Coulomb 
potential. Both types of particles are explicitly treated 
as paths and all correlation effects are included, which 
makes PIMC one of the most accurate finite-temperature 
quantum simulation methods available. The only non- 
controlled approximation required is the fixed node ap- 
proximation that is introduced to treat the fermion sign 
problem, which arises from the explicit treatment of elec- 
trons. The fermion nodes are taken from a thermal trial 
density matrix, for which we extended the variational 
density matrix approach |l7| to helium. 

We complement our PIMC EOS with DFT-MD data, 
since DFT is much more efficient at low temperatures 
because it is in principle a ground state electronic struc- 
ture method. The DFT-MD trajectories were obtained 
with Born-Oppcnheimer MD where the electrons were as- 
sumed to be in the instantaneous ground state. We used 
the CPMD code llSl using the PBE generalized gradi- 
ent approximation [19| | with N—64 atoms, a time step of 
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FIG. 1: The principal shock Hugoniot curves computed using 
DFT-MD with (solid line) and without (blue dashed line) the 
consideration of excited electronic states are compared with 
gas gun shock wave experiments pH (open red symbols) . The 
dashed purple line shows the computed reshock curve. The o 
on it indicates the impedance match condition for the reshock 
experiment pH (solid red symbol). 

0.77 fs, Troullier-Martin pseudopotentials, and a 100 Ry 
cut-off for the plane wave expansion of the Kohn Sham 
orbitals, combined with T point sampling of the Brillioun 
zone. A finite size study showed that the Hugoniot re- 
sults are well converged with N=6A. 

The only available shock data for fluid helium were 
obtained with gas gun experiments by Nellis et al. |l6|. 
The comparison shown in Fig. ^ shows excellent agree- 
ment between experimental data and DFT-MD simula- 
tion results for the principal and the reshock Hugoniot 
curves. 

The DFT-MD results without electronic excitations 
are closely tracked by data from classical Monte Carlo 
(CMC) simulations using the Aziz pair potential [2p| . 
This potential was derived to describe the interaction of 
two isolated helium atoms. With rising shock pressure, 
DFT-MD and CMC curves in Fig. show a gradual in- 
crease in compression towards 4-fold the initial density, 
which represents the high pressure limit, where the ki- 
netic energy dominates over the potential energy and the 
system behaves approximately like non-interacting par- 
ticles. 

Conversely, the results from PIMC calculations predict 
much higher compression ratios, reaching a maximum of 
5.24±0.04 for P = 360 GPa and T = 150 000 K. This 
increase in compression is due to electronic excitations 
that are present in the hot dense fluid at thermodynamic 
equilibrium. Helium differs in this regard from shocked 
deuterium, for which PIMC calculation predicted a max- 
imum compression ratio of only 4.3(1) |8|]. 

In our PIMC simulation program, the pair density ma- 




FIG. 2: The principal shock Hugoniot curve of helium is 
shown as a function of compression ratio. PIMC and DFT- 
MD results were combined to span a large range of temper- 
ature, as the circles in the curve indicate. Effects of thermal 
electronic excitations (purple solid line) were used to correct 
the ground state DFT-MD calculations (blue dashed line), 
which resulted in a strong increase in compressibility. Re- 
sults a from classical simulation using the Aziz pair potential 
and from a non-interacting plasma model are included. 

trix technique [2lL E3 | is employed to treat the Coulomb 
interactions. The accurate treatment of the many-body 
correlations and the application of the nodal restriction 
for the paths requires the use of a small step, r = /3/M, 
to discretize the paths in imaginary time /3 = l/fceT into 
M steps. Using a time step between r _1 = 2 x 10 6 and 
16 x 10 6 K depending on density, allowed a reduction of 
the remaining time step error in the observables below 
the size of MC error bars. 

PIMC simulations with N e — 32 electrons and iV=16 
nuclei were performed on a grid of 12 densities and 5 tem- 
peratures ranging from 3.5 > r s > TO and 10 6 > T > 
61250 K. The resulting EOS was interpolated to obtain 
the Hugoniot curves. A finite-size extrapolation with up 
to N—57 nuclei was done for four points in T-p space, 
T=125 000 and 10 6 K combined with r s =1.86 and 1.25. 
Above 100 GPa, the Hugoniot curve is insensitive to fi- 
nite size corrections because of the high temperatures 
and due to the partial cancellation of pressure and inter- 
nal energy corrections in Eq. ^ After propagating the 
correction and the uncertainty of the finite size extrapo- 
lation, we obtained 5.24±0.04 as maximum compression 
ratio of fluid helium, which actually brackets the original 
value of 5.25 obtained for N—1Q. 

The predicted increase in compression beyond 4-fold 
the initial density in helium can be understood by invok- 
ing a very simple free energy model, which assumes he- 
lium to be composed of different non-interacting species, 
He, He + , and He 2+ , to represent the various ionization 
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stages of helium as well as free electrons. The resulting 
Hugoniot curve shows a similar increase in the compres- 
sion ratio demonstrating that the increase beyond 4-fold 
is due to electronic excitations leading to free electrons. 

To further verify this hypothesis, we corrected the DFT 
EOS for finite-temperature electronic effects. For a num- 
ber of snapshots along the MD trajectory we therm ally 
populated the instantaneous excited electronic states |23j 
using the Mermin functional with up to 7 additional or- 
bitals per atom. For temperatures above 15 000 K, the 
resulting corrections to the internal energy and pressure 
leads to a substantial increase in shock compressibility. 
This increase is primarily caused by a rise in the inter- 
nal energy due to thermal population of excited elec- 
tronic states. Fluid helium maintains a wide excitation 
gap ranging between 5 and 15 eV for T < 80 000 K and 
2.6 > r s > 1.0. On the principal Hugoniot, electronic 
excitations occur above 20 GPa, which explains why the 
gas gun experiments have not reached the regime of elec- 
tronic excitations. Even a gas gun reshock experiment 
would be insufficient because the final temperatures re- 
main relatively low (Fig. ^) , and facilities that can gen- 
erate faster shock waves are needed instead. 

The maximum shock pressures that can be reached at 
a particular experimental facility depend on the power of 
the drive but also on the impedance of the sample mate- 
rial. For the same shock drive, the final shock pressure 
is, to a first approximation, proportional to the initial 
density of the sample material. The maximum pressure 
reported for deuterium Nova laser shocks [10, 340 GPa, 
and the highest achievable on the Z machine 175 
GPa, would translate to approximately 246 GPa, and 
127 GPa [l^l in helium. Consequently, with both facili- 
ties one would be able reach the regime of the predicted 
5-fold compression and probe the effect of electronic ex- 
citations. 

The comparison in Fig. |31 shows the discussed increase 
in compressibility beyond 4-fold for helium, while our 
results for deuterium ||, |2(| show hardly any, despite 
the similarity of the two fluids. Deuterium molecules 
and helium atoms have the same mass and provide two 
mechanisms to absorb shock energy, which in principle 
can lead to shock compression ratios substantially larger 
than 4. The helium atom has two ionization stages with 
energies of 24.6 and 54.4 eV. Deuterium molecules disso- 
ciate with 4.5 eV energy, and the ionization of resulting 
atoms requires 13.6 eV. However, the explanation for the 
different shock behavior of helium and deuterium is not 
a consequence of single particle properties but is a re- 
sult of different degrees of particle interaction. Fig. 0] 
shows the Hugoniot function, H, for both materials at 
5-fold compression. For helium, first principles calcu- 
lation that include the interaction effects, as well as the 
non-interacting plasma model, predict that H changes its 
sign, implying that helium is more than 5-fold compress- 
ible. While the non-interacting plasma model predicts 
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FIG. 3: The principal shock Hugoniot curves for deu- 
terium [2^| and helium are shown for samples that were pre- 
compressed to different initial densities. For both materials, 
the precompression reduces the maximum compression ratio 
p/po that can be reached. For helium, the indicated initial 
densities (a)-(d) correspond to the initial pressures of 7.1 kbar, 
1.8 kbar, 188 bar |25|. and 1 bar. The DFT-MD and PIMC 
data shown in Fig. |5|were interpolated for helium. 

the same behavior also for deuterium, the PIMC simu- 
lations show that the H function is strictly negative be- 
cause the pressure is substantially higher than suggested 
by the non-interacting plasma model. 

The interaction effects at 5-fold compression are 
stronger in deuterium than in helium because in helium, 
the initial density is lower and, more importantly, the 
electronic orbitals are much more localized. The initial 
density, po = 0.1235 gem -3 corresponds to a Wigner- 
Seitz radius of r s = 3.51 [f 7r(r s a ) 3 = V/N B ], and 95% of 
the electron charge is localized in 10% of the total volume. 
The initial density of fluid deuterium, po = 0.171 gem 
corresponds to r s — 3.16, and 95% of the charge is local- 
ized in 40% of the volume. With increasing compression, 
the orbitals in deuterium start to overlap sooner and the 
resulting polarization and exchange interactions will in- 
crease the pressure, which causes deuterium to appear 
less compressible in shock experiments. 

In Fig. El we use our first principles calculations to pre- 
dict how the static precompression will affect the Hugo- 
niot curves of helium and deuterium |24| . For helium, 
a precompression to 4-fold the cryogenic liquid density 
at 1 bar, po = 4 / 0Q lbar \ increases the maximum density 
on the Hugoniot curve to p = 17.5 p^ hal \ However, the 
maximum compression ratio for the dynamic compres- 
sion, r\ = p/po, decreases in both materials with increas- 
ing precompression. The increase in the initial density 
leads to stronger interactions under shock conditions. In 
general the shock compression ratio is determined by the 
relative importance of excitations of internal degrees of 
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FIG. 4: The Hugoniot functions H from Eq.^for helium and 
deuterium derived from first principles calculation (solid lines) 
and from the non-interacting plasma model (dashed line) is 
shown at 5-fold compression. Deuterium is less compressible 
than helium because the interaction between the particles is 
much stronger, which leads to a substantial increase in pres- 
sure beyond the corresponding value of the non-interacting 
system including ionization (-Pid), illustrated by the pressure 
difference in the lower graph. 



freedom that increase the internal energy and interaction 
effects that increase the pressure. 

The point of maximum compression, rj max , along the 
Hugoniot curve is reached when the Griineisen parame- 
ter, 7 = V satisfies 7 = A reduction in the 
compression ratio with precompression can be expressed 
by d ^f x > 0. Using Eq. ^ an d assuming Eq does not 



dV 

change with the precompression, this condition can be 
expressed in terms of the simple thermodynamic condi- 



tion, 



OP 

dp 



> 1. We have computed the isoenergetic 



compressibility, 



and verified that this condition is 



satisfied for hydrogen and helium by using our first prin- 
ciples results as well as with the non-interacting plasma 
model. This confirms that the precompression reduces 
the maximum compression ratio in hydrogen and helium, 
and one can postulate that this might also be true for 
other simple fluids. 

The computed EOS shows that helium and hydrogen 
behave differently at high pressure, which has important 
implications for the interior structure of solar and extra- 
solar giant planets where it is has been predicted that the 
two fluid phases become immiscible at high pressures [27j . 
Our calculations predict that helium is 5.24(4)-fold com- 
pressible under shock conditions, which distinguishes it 
from fluid deuterium, for which firstprinciples calcula- 
tions and recent experiments [3J, U, |5j, |(| predicted a max- 
imum compression ratio close to 4.3. We suggest that all 



deuterium experiments be repeated with helium in order 
to validate the different shock compression techniques. 
This validation is important before the EOS can be used 
for wide range of applications including the modeling of 
giant planets and to draw conclusions about their evolu- 
tion. 
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